#ifndef SPHEREWRAPPER_H_
#define SPHEREWRAPPER_H_

#include <vector>
#include <math.h>

#include "olb3D.h"
#include <ibm/ibm3D.h>

class IntegrationPoint
{
public:
	IntegrationPoint(double _x=0.0, double _y=0.0, double _z=0.0, double integrationArea=0.0, double integrationLength=0.0)
		: x(_x), y(_y), z(_z), intArea(integrationArea), intLen(integrationLength) {}
	
	double x,y,z;
	
	void SetIntegrationArea(double intArea) { this->intArea = intArea; }
	double GetIntegrationArea() const { return this->intArea;}	

	void SetIntegrationLength(double intLength) { this->intLen = intLength; }
	double GetIntegrationLength() const { return this->intLen;}	

	void        SetCurrentInteractionForce(NJRvector3d f) { this->force = f; }
	NJRvector3d GetCurrentInteractionForce() const { return this->force; }
	
	double GetIntegrationVolume() const { return this->intArea * this->intLen;}	

	void SetNormal(NJRvector3d normal) { this->normal = normal; }
	NJRvector3d GetNormal() const { return this->normal;}	

private:
	double intArea, intLen;
	NJRvector3d normal, force;
	int doID;
};

class SphereWrapper
{
public:
	SphereWrapper(int strip=10, double thick=1.0); 
	virtual ~SphereWrapper() {}
	void Reset();
	void WrapSphere(double x, double y, double z, double r, double rotateMatrix[3][3]);
	std::vector<IntegrationPoint*>& GetIntegrationPoints() { return this->currIntPts; }
	
private:
	void RecreateIntegrationPoints(int strip, double thick);
	std::vector<IntegrationPoint*> currIntPts, initIntPts;
	double center_x, center_y, center_z;
};

SphereWrapper::SphereWrapper(int strip, double thick)
{ 
	this->RecreateIntegrationPoints(strip, thick); 
}

void SphereWrapper::RecreateIntegrationPoints(int strip, double thick)
{
	this->initIntPts.clear();
	this->currIntPts.clear();
	double radius = 1.0;
	double bandWidth = M_PI / (double)strip;
	double rr = radius*radius;
	for(int pt=0; pt<strip+1; pt++)
	{
		if(pt == 0)
		{
			double x = 0.0; double y = 0.0; double z = -radius;
			double theata1 = -0.5 * M_PI;
			double theata2 = theata1 + 0.5*	bandWidth;
			double area = 2.0 * M_PI * rr * (sin(theata2) - sin(theata1));
			this->initIntPts.push_back(new IntegrationPoint(x,y,z,area,thick));		
			this->currIntPts.push_back(new IntegrationPoint(x,y,z,area,thick));		
		}
		else if(pt == strip)
		{
			double x = 0.0; double y = 0.0; double z = radius;
			double theata2 = 0.5 * M_PI;
			double theata1 = theata2 - 0.5*	bandWidth;;
			double area = 2.0 * M_PI * rr * (sin(theata2) - sin(theata1));
			this->initIntPts.push_back(new IntegrationPoint(x,y,z,area,thick));		
			this->currIntPts.push_back(new IntegrationPoint(x,y,z,area,thick));		
		}
		else
		{
			double theata = -0.5 * M_PI + bandWidth * pt;
			double z = radius * sin(theata);
			double theata1 = theata - 0.5 * bandWidth;
			double theata2 = theata + 0.5 * bandWidth;
			double deltaArea = 2.0 * M_PI * rr * (sin(theata2) - sin(theata1)) / (double)(strip);
			double deltaPhi = 2.0 * M_PI / (double)(strip);
			for(int pt2=0; pt2<strip; pt2++)
			{
				double phi = deltaPhi * pt2;
				double x = radius * cos(phi); double y = radius * sin(phi);
				this->initIntPts.push_back(new IntegrationPoint(x,y,z,deltaArea,thick));		
				this->currIntPts.push_back(new IntegrationPoint(x,y,z,deltaArea,thick));		
			}
		}
	}
	this->Reset();
}

void SphereWrapper::Reset()
{
	this->center_x = 0.0;
	this->center_y = 0.0;
	this->center_z = 0.0;
	for(unsigned int i=0; i < this->currIntPts.size(); i++)
	{
		this->currIntPts[i]->x = this->initIntPts[i]->x;
		this->currIntPts[i]->y = this->initIntPts[i]->y;
		this->currIntPts[i]->z = this->initIntPts[i]->z;
		this->currIntPts[i]->SetIntegrationArea(this->initIntPts[i]->GetIntegrationArea());
		this->currIntPts[i]->SetIntegrationLength(this->initIntPts[i]->GetIntegrationLength());
		NJRvector3d vec(this->initIntPts[i]->x, this->initIntPts[i]->y, this->initIntPts[i]->z);
		NJRvector3d normal = vec * (1./vec.length());
		this->currIntPts[i]->SetNormal(normal);
	}
}

void SphereWrapper::WrapSphere(double x, double y, double z, double r, double rotMat[3][3])
{
	this->Reset();
	//set Radius
	for(unsigned int i=0; i < this->currIntPts.size(); i++)
	{
		this->currIntPts[i]->x *= r;
		this->currIntPts[i]->y *= r;
		this->currIntPts[i]->z *= r;
		this->currIntPts[i]->SetIntegrationArea(r*r*this->initIntPts[i]->GetIntegrationArea());
	}

	//rotate
	if(rotMat)
	{
		for(unsigned int i=0; i < this->currIntPts.size(); i++)
		{
			double x = this->currIntPts[i]->x;
			double y = this->currIntPts[i]->y;
			double z = this->currIntPts[i]->z;
			double x1 = rotMat[0][0] * x + rotMat[0][1] * y + rotMat[0][2] * z;  
			double y1 = rotMat[1][0] * x + rotMat[1][1] * y + rotMat[1][2] * z;  
			double z1 = rotMat[2][0] * x + rotMat[2][1] * y + rotMat[2][2] * z;  	
			this->currIntPts[i]->x = x1;
			this->currIntPts[i]->y = y1;
			this->currIntPts[i]->z = z1;
			NJRvector3d vec(this->currIntPts[i]->x - this->center_x, this->currIntPts[i]->y - this->center_y,  this->currIntPts[i]->z - this->center_z);
			NJRvector3d normal = vec * (1./vec.length());
			this->currIntPts[i]->SetNormal(normal);
		}
	}

	//translate
	for(unsigned int i=0; i < this->currIntPts.size(); i++)
	{
		this->currIntPts[i]->x += x;
		this->currIntPts[i]->y += y;
		this->currIntPts[i]->z += z;
	}
}

#endif /*SPHEREWRAPPER_H_*/
